#!/bin/sh
#
# Run FDR analysis
# 
# by Marcelo G Mattar (08/19/2012)
# mattar@sas.upenn.edu

STATSDIR=$1
COPE=$2
FDRRATE=$3

echo
echo "Working on $STATSDIR; Contrast #$COPE"

# Convert T-values to P-values using the cope and variance cope images
ttologp -logpout $STATSDIR/logp${COPE} $STATSDIR/varcope${COPE} $STATSDIR/cope${COPE} `cat $STATSDIR/dof`

# Caclulate the p-value image
fslmaths $STATSDIR/logp${COPE} -exp $STATSDIR/p${COPE}

# Calculate the FDR from the p-value image
fdr -i $STATSDIR/p${COPE} -m $STATSDIR/../mask -q $FDRRATE > $STATSDIR/fdr${COPE}

# Get the FDR threshold
echo "FDR (rate) is" $FDRRATE
FDRTHRESH=`awk 'NR==2' $STATSDIR/fdr${COPE}`
echo "FDR p-threshold is" $FDRTHRESH
FDRTHRESH=`awk -v a=$FDRTHRESH -v b="1" 'BEGIN{print (b - a)}'`
echo "1 minus FDR threshold is" $FDRTHRESH

# Calculate a p-value image using the FDR threshold 
fslmaths $STATSDIR/p${COPE} -mul -1 -add 1 -thr $FDRTHRESH -mas $STATSDIR/../mask $STATSDIR/thresh_fdr_1_minus_p${COPE}

# Calculate how many voxels survived the applied threshold
NUMVOXS=`fslstats $STATSDIR/thresh_fdr_1_minus_p${COPE} -V`
NUMVOXS=`echo $NUMVOXS | cut -d' ' -f1`

echo
echo "$NUMVOXS voxels survived the FDR threshold"

# Clean up unneccesary files
#rm -f $STATSDIR/logp${COPE}.nii.gz $STATSDIR/p${COPE}.nii.gz fdr${COPE}